function handles = analyzeEPArray_detectAPs(handles)
% analyzeEPArray_detectAPs     
%
% initial, detection-only part of original function
%
% gs june 2006
% ------------------------------------------------------------------------------

handles.data.analysis.cellAttachedCheck = 1; 
% Note: this is basically a flag, indicating that we're in EP analysis mode

dataArray = handles.data.map.mapActive.bsArray; % dataArray
% TODO: eventually should probably make this conditional ...

[rows, numTraces] = size(dataArray); % must be 2D data array; col num = trace num
sr = handles.data.acq.sampleRate;
mapPattern = handles.data.map.mapActive.mapPattern;
traceData = [];

% TIME vector
totalTime = (rows-1)/sr; % in sec, not msec
xTimeAxis = linspace(0, totalTime, rows)';
handles.data.map.mapActive.xTimeAxis = xTimeAxis;
minX = 0; 
maxX = max(xTimeAxis);

% BASELINE settings, mean, sd, baseline-subtracted array
baselineStartIndex = round(handles.data.analysis.baselineStart*sr)+1;
baselineEndIndex = round(handles.data.analysis.baselineEnd*sr)+1;
baselineMeans = mean(dataArray(baselineStartIndex:baselineEndIndex,:));
handles.data.map.mapActive.baselineMeans = baselineMeans';
baselineSDs = std(dataArray(baselineStartIndex:baselineEndIndex,:));
handles.data.map.mapActive.baselineSDs = baselineSDs';
bsArray = handles.data.map.mapActive.bsArray;

% NOISE of baseline
%       b0 = 0-50 msec, baseline for baseline window
%       r0 = 50-100 msec, baseline window
%       b1 = 50-100 msec, baseline for synaptic window
%       r1 = 110-160 msec, synaptic window
% % mean 
% b0 = mean(state.uncaging.map.bsArray(1:500, :));
% b1 = mean(state.uncaging.map.bsArray(500:1000, :));
% r0 = mean(state.uncaging.map.bsArray(500:1000, :));
% r1 = mean(state.uncaging.map.bsArray(1101:1601, :));
% peak 
b0 = min(handles.data.map.mapActive.bsArray(1:500, :));
b1 = min(handles.data.map.mapActive.bsArray(500:1000, :));
r0 = min(handles.data.map.mapActive.bsArray(500:1000, :));
r1 = min(handles.data.map.mapActive.bsArray(1101:1601, :));
% TODO: make the timing set-able from the gui, or automatically based on the baseline window params
handles.data.map.mapActive.ampHistoBaselineData = -(r0-b0);
handles.data.map.mapActive.histoBaseSDs = std(handles.data.map.mapActive.ampHistoBaselineData);
handles.data.map.mapActive.ampHistoSynapticData = -(r1-b1);

% STIMULUS timing
stimOn = handles.data.analysis.stimOnEP;
stimOnInd = round(stimOn*sr+1);
handles.data.analysis.stimOnInd = stimOnInd;
% stimOff = stimOn + handles.data.analysis.shutterOpen;                
% stimOffInd = round(stimOff*sr+1);

% WINDOWS *******************************************************************
directWindowStart = handles.data.analysis.responseStartEP;
directStartInd = round(directWindowStart * sr + 1);
directWindow = handles.data.analysis.responseDurEP;
directEndInd = directStartInd + round(directWindow * sr);
synStart = directWindow;
synStartInd = directEndInd;
% synapticWindow = handles.data.analysis.synapticWindow;
% synEndInd = directStartInd + round(synapticWindow * sr);
synEndInd = directEndInd;

% THRESHOLDS ****************************************************************
% NB: these are vectors
handles.data.map.mapActive.apThreshold = handles.data.analysis.apThreshold * ones(1, numTraces)';

% distance array
[R, C] = size(handles.data.map.mapActive.mapPattern);
spacingXY = [handles.data.map.mapActive.mapSpacing handles.data.map.mapActive.mapSpacing];
try
    somaXYZ = [handles.data.acq.somaXnew handles.data.acq.somaYnew handles.data.acq.somaZ];
catch
    somaXYZ = [handles.data.acq.somaX handles.data.acq.somaY handles.data.acq.somaZ];
    disp('analyzeEPArray: using original (untransformed) soma X and Y coordinates; check that this is accurate');
end
% handles.data.map.mapActive.distFromSoma = calcDistFromSoma(spacingXY, R, C, somaXYZ);
% handles.data.map.mapActive.distFromSoma = reshape(handles.data.map.mapActive.distFromSoma, [], 1);

% EVENTS ****************************************************************
% APs
% if handles.data.analysis.currentClamp | handles.data.analysis.cellAttachedCheck % shouldn't need this; GS20060530
    if ~handles.data.analysis.cellAttachedCheck
        result = arrayEventFinder(dataArray(stimOnInd:synEndInd,:), handles.data.map.mapActive.apThreshold, 'positive');
    else
        result = arrayEventFinder(bsArray(stimOnInd:synEndInd,:), handles.data.map.mapActive.apThreshold, 'positive');
    end
    events = zeros(numTraces,1);
    if ~isempty(result)
        tracesWithEvents = result(:,1);
        timeIndicesOfEvents = result(:,2);
        onsetLatEvents = result(:,2)/sr;
        events(tracesWithEvents) = onsetLatEvents;
    end
    handles.data.map.mapActive.latAPonset = events;
    tracesWithAPs=find(handles.data.map.mapActive.latAPonset);
%     disp(' ');
%     disp(['     Traces with APs:     ' num2str(tracesWithAPs')]);
%     disp(' ');
    if ~isempty(result) % i.e., if there were APs
        nAPmax = 8; % max number of APs in latency arrays -- 
        % TODO: make this GUI-settable
        handles.data.map.mapActive.apNumber = zeros(numTraces,1);
        handles.data.map.mapActive.apDelay = zeros(numTraces,nAPmax); 
        for n = 1:length(tracesWithAPs) % process one trace at a time
            nAP=0; AP=0; peak=0; X=0; apDelay = [];

            if ~handles.data.analysis.cellAttachedCheck
                traceData = dataArray(stimOnInd:synEndInd, tracesWithAPs(n));
            else
                traceData = bsArray(stimOnInd:synEndInd, tracesWithAPs(n));
            end
            
            for i = 2:length(traceData)-1
                if traceData(i) >= handles.data.analysis.apThreshold
                    AP = 1;
                    if traceData(i-1) <= handles.data.analysis.apThreshold
                        spikeStart = i;
                        nAP = nAP+1;
                    end
                    if traceData(i+1) < handles.data.analysis.apThreshold & exist('spikeStart')
                        spikeEnd = i;
                        [peak(nAP), X(nAP)] = max(traceData(spikeStart:spikeEnd));
                        X(nAP) = (X(nAP) + spikeStart-2)/sr; % points -> time
                        apDelay(nAP) = X(nAP);
                    end
                end
            end
            
            handles.data.map.mapActive.apNumber(tracesWithAPs(n)) = nAP;
            handles.data.map.mapActive.apDelay(tracesWithAPs(n), 1:length(apDelay)) = apDelay;
        end
        
        handles.data.map.mapActive.apBinary = handles.data.map.mapActive.apNumber > 0;
        handles.data.map.mapActive.apDelayArray = handles.data.map.mapActive.apDelay(tracesWithAPs,:);
    else
        handles.data.map.mapActive.apNumber = 0;
        handles.data.map.mapActive.apBinary = 0;
        handles.data.map.mapActive.apDelay = [];
        handles.data.map.mapActive.apDelayArray = [];
    end
        
% end % shouldn't need this; GS20060530

% % spit out some AP results
% if handles.data.analysis.cellAttachedCheck
%     disp('% EXCITATION PROFILE ANALYSIS:');
%     disp(['map.tracesWithSpikes = [' num2str(tracesWithAPs') '];']);
%     disp(['map.numberOfSpikes = [' num2str(handles.data.map.mapActive.apNumber(handles.data.map.mapActive.apNumber>0)') '];']); 
%     disp(['map.totalNumberOfSites = ' num2str(length(tracesWithAPs')) ';']);
%     disp(['map.totalNumberOfSpikes = ' num2str(sum(handles.data.map.mapActive.apNumber(handles.data.map.mapActive.apNumber>0)')) ';']);
%     if length(tracesWithAPs)>0
%         spikesPerSite = num2str(round(100*sum(handles.data.map.mapActive.apNumber(handles.data.map.mapActive.apNumber>0)')/length(tracesWithAPs'))/100);
%         meanDist = num2str(round(mean(handles.data.map.mapActive.distFromSoma(handles.data.map.mapActive.apBinary'))));
%     else
%         spikesPerSite = 'NaN';
%         meanDist = 'NaN';
%     end
%     disp(['map.spikesPerSite = ' spikesPerSite ';']);
%     disp(['map.meanDistanceFromSoma = ' meanDist ';']);
%     disp('map.spikeLatencies = [');
%     disp(handles.data.map.mapActive.apDelay(tracesWithAPs,:));
%     disp('];');    
%     handles.data.analysis.currentClamp = 1; % may want to improve this kluge eventually
% end

% MAP display ******************************************************************

handles = makeMapDataArray(handles);
disp('% EXCITATION PROFILE ANALYSIS:');
disp('map.apBinary = [');
disp(handles.data.map.mapActive.mapDataArray(:,:,8));
disp('];');
disp('map.apNumber = [');
disp(handles.data.map.mapActive.mapDataArray(:,:,6));
disp('];');
disp('map.apDelay = [');
disp(handles.data.map.mapActive.mapDataArray(:,:,7));
disp('];');
disp('map.spikeLatencies = [');
disp(handles.data.map.mapActive.apDelayArray);
disp('];');    

% gs 20060616
% % % maps = handles.data.map.mapActive.mapDataArray(:,:,6);
% % % spacingXY = 
% % % somaXYZ = 
% % figFlag = 1;
% % handles = excitationProfileAnalysis(handles, spacingXY, somaXYZ, figFlag);
% % 
% % try
% %     handles = array2DplotEP(handles);
% % catch
% % end


% GS20060531 - commented; no longer needed (& causing error)
% % % ---------------------------------------------------------------------
% % function out = calcDistFromSoma(spacingXY, R, C, somaXYZ)
% % 
% % % arrays and vectors of distances
% % colVals = [-((R-1)/2)*spacingXY(1):spacingXY(1):((R-1)/2)*spacingXY(1)]';
% % colVals = colVals - somaXYZ(1);
% % colArray = repmat(colVals, 1, C);
% % 
% % rowVals = [-((C-1)/2)*spacingXY(2):spacingXY(2):((C-1)/2)*spacingXY(2)];
% % rowVals = rowVals - somaXYZ(2);
% % rowArray = repmat(rowVals, R, 1);
% % 
% % planeArray = somaXYZ(3) * ones(R,C);
% % 
% % distArray = sqrt(rowArray.^2 + colArray.^2 + planeArray.^2);
% % % dists = round(distArray), return
% % [r,c] = size(distArray);
% % distances = reshape(distArray, r*c, []);
% % distances2 = [-distances; distances];
% % 
% % out = distances;